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Abstract 

The knowledge of transitions between regular, laminar or chaotic behavior is essential to 
understand the underlying mechanisms behind complex systems. While several linear ap- 
proaches are often insufficient to describe such processes, there are several nonlinear meth- 
ods which however require rather long time observations. To overcome these difficulties, we 
propose measures of complexity based on vertical structures in recurrence plots and apply 
them to the logistic map as well as to heart rate variability data. For the logistic map these 
measures enable us not only to detect transitions between chaotic and periodic states, but 
also to identify laminar states, i. e. chaos-chaos transitions. The traditional recurrence quan- 
tification analysis fails to detect the latter transitions. Applying our new measures to the 
heart rate variability data, we are able to detect and quantify the laminar phases before a 
life-threatening cardiac arrhythmia occurs thereby facilitating a prediction of such an event. 
Our findings could be of importance for the therapy of malignant cardiac arrhythmias. 
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I. INTRODUCTION 



Numerous scientific disciplines, such as astrophysics, biology or geosciences, use 
data analysis techniques to get an insight into the complex processes observed in 
nature ^ 0] which show generally a nonstationary and complex behavior. As 
these complex systems are characterized by different transitions between regular, 
laminar and chaotic behaviors, the knowledge of these transitions is necessary for 
understanding the process. However, observational data of these systems are typi- 
cally rather short. Linear approaches of time series analysis are often not sufficient 

^] and most of the nonlinear techniques (cf [j^, |7p), such as fractal dimensions or 
Lyapunov exponents []7|, suffer from the curse of dimensionality and require 

rather long data series. The uncritical application of these methods, especially to 
natural data, can therefore be very dangerous and it often leads to serious pitfalls. 

To overcome these difficulties other measures of complexity have been proposed, 
such as the Renyi entropies, the effective measure complexity, the ^-complexity or 
the renormalized entropy \^^, |T^]. They are mostly based on symbolic dynamics and 
are efficient quantities for characterizing measurements of natural systems, such as 
in cardiology []T^, [T^, 0, cognitive psychology [ [lEQ or astrophysics []T7|, |T^, [T^. In 



this paper we focus on another type of measures of complexity, which is based on the 
method of recurrence plots (RP). This approach has been introduced for the analysis 
of nonstationary and rather short data series [ pD| , PT| , |2^ ]. Moreover, a quantitative 
analysis of recurrence plots has been proposed to detect typical transitions (e. g. bi- 
furcation points) occurring in complex systems ] ^3\ , |2^ , ^SQ . However, the quantities 
introduced so far are not able to detect more complex transitions, especially chaos- 
chaos transitions, which are also typical in nonlinear dynamical systems. Therefore 
in this paper we introduce measures of complexity based on recurrence plots which 
allow us to identify laminar states and their transitions to regular as well as other 
chaotic regimes in complex systems. These measures make the investigation of in- 
termittency of processes possible, even if they are only represented by short and non- 
stationary data series. 

The paper is organized as follows: After a short review of the technique of re- 
currence plots and some measures we introduce new measures of complexity based 
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on recurrence plots. After that we apply the new approach to the logistic equation 
and demonstrate the ability to detect chaos-chaos transitions. Finally, we apply this 
technique to heart rate variability data [pEp. We demonstrate that by applying our 
new proposed methods we are able to detect laminar phases before the onset of a 
life-threatening cardiac arrhythmia. 



II. RECURRENCE PLOTS AND THEIR QUANTIFICATION 

The method of recurrence plots (RP) was firstly introduced to visualize the time 
dependent behavior of the dynamics of systems, which can be pictured as a trajectory 
Xi e TZ^ (i = 1, . . . , N) in the n-dimensional phase space [piQ. It represents the recur- 
rence of the phase space trajectory to a certain state, which is a fundamental property 
of deterministic dynamical systems []27|, The main step of this visualization is the 



calculation of the N x A^-matrix 



Rij:=e{ei-\\xi-Xj\\), i,j = l...N, (1) 

where ei is a cut-off distance, || ■ || a norm (e.g. the Euclidean norm) and 6(x) 
the Heaviside function. The phase space vectors for one-dimensional time series Ui 
from observations can be reconstructed by using the Taken's time delay method Xi = 
{ui, Ui+r, ■ ■ ■ , Mi+(m-i) r) The dimeusiou m can be estimated with the method of false 
nearest neighbours (theoretically, m = 2n + 1) [0, The cut-off distance Si defines 
a sphere centered at Xj. If £, falls within this sphere, the state will be close to Xi and 
thus Rj j = 1. These Si can be either constant for all xi [ |22[ j or they can vary in such a 



way, that the sphere contains a predefined number of close states []2ip. In this paper 
a fixed Si and the Euclidean norm are used, resulting in a symmetric RP. The binary 
values in Rij can be simply visualized by a matrix plot with the colors black (1) and 
white (0). 

The recurrence plot exhibits characteristic large-scale and small-scale patterns 



which are caused by typical dynamical behavior [^, §4p, e. g. diagonals (similar local 
evolution of different parts of the trajectory) or horizontal and vertical black lines 
(state does not change for some time). 
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Zbilut and Webber have recently developed the recurrence quantification analy- 
sis (RQA) to quantify an RP ^]. They define measures using the recur- 
rence point density and the diagonal structures in the recurrence plot, the recurrence 
rate, the determinism, the maximal length of diagonal structures, the entropy and 
the trend. A computation of these measures in small windows moving along the main 
diagonal of the RP yields the time dependent behavior of these variables and, thus, 
makes the identification of transitions in the time series possible [^]. 

The RQA measures are mostly based on the distribution of the length of the diago- 
nal structures in the RP. Additional information about further geometrical structures 
such as vertical and horizontal elements are not included. Gao has therefore recently 
introduced a recurrence time statistics, which corresponds to vertical structures in 
an RP I P^I , In the following, we extend this view on the vertical structures and 



define measures of complexity based on the distribution of the vertical line length. 
Since we are using symmetric RPs here, we will only consider the vertical structures. 

III. MEASURES OF COMPLEXITY 

We consider a point Xj of the trajectory and the set of its associated recurrence 
points Si := {xk : Ri,fc = 1 ; k e [1 . . . N — 1]}. Denote a subset of these recurrence 
points Si := {xi G Si : (R^,; ■ Ri,«+i) + (Ri,; • Rij-i) > ; / G [1 . . . iV], Ri,o = Ri,7V+i := 0} 
which contains the recurrence points forming the vertical structures in the RP at 
column i. In continuous time systems with high time resolution and with a not too 
small threshold e, a large part of this set Si usually corresponds to the sojourn points 
described in [ ^9] , Although sojourn points do not occur in maps, the subset Si 
is not necessarily empty. Next, we determine the length v of all connected subsets 
{xj ^ Si] Xj+i . . . Xj+i, G Si] ^ Si} in Si. Piiv) = [vi] I = 1,2, ... L} denotes the set 

of all occurring subset lengths in Sj and from IJ^^^ Pi{v) we determine the distribution 
of the vertical line lengths P{v) in the entire RP. 



Analogous to the definition of the determinism ^Tp, we compute the ratio be- 



tween the recurrence points forming the vertical structures and the entire set of re- 
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currence points 

A := ^"="""" (2) 

and call it laminarity A. The computation of A is realized for v which exceeds a 
minimal length Vmin- For maps we use Vmin = 2. A is the measure of the amount of 
vertical structures in the whole RP and represents the occurrence of laminar states in 
the system, without, however, describing the length of these laminar phases. It will 
decrease if the RP consists of more single recurrence points than vertical structures. 
Next, we define the average length of vertical structures 

vPiv) 

what we call trapping time T. The computation also uses the minimal length Vmin as 
in A. The measure T contains information about the amount and the length of the 
vertical structures in the RP. 

Finally, we use the maximal length of the vertical structures in the RP 

Vmax = max {{vi ; / = 1, 2, . . . L}) (4) 

as a measure, which is the analogue to the standard RQA measure Lmax D^]- 

Although the distribution of the diagonal line lengths also contains information 
about the vertical line lengths, the two distributions are significantly different. In 
order to compare the measures proposed with the standard RQA measures, we apply 
them to the logistic map. 



IV. APPLICATION TO THE LOGISTIC MAP 



In order to investigate the potentials of A, T and Vmax, we firstly analyze the logistic 
map 

Xn+l = Oa^n (1 ~ Xn) , (5) 

especially the interesting range of the control parameter a e [3.5, 4] with a step width 
of Aa = 0.0005. Starting with the idea of TruUa et al. [ ^BQ to look for vertical struc- 
tures, we are especially interested in finding the laminar states in chaos-chaos tran- 
sitions. Therefore we generate for each control parameter a a separate time series. In 
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the analyzed range of a G [3.5, 4] various regimes and transitions between them occur, 
e.g. accumulation points, periodic and chaotic states, band merging points, period 
doublings, inner and outer crisis 



Logistic map 
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FIG. 1: (A) Bifurcation diagram of the logistic map. (B) Low ordered supertrack functions 
Si{a) (i = 1 . . . 10) and the fixed point of the logistic map 1 — 1/a (dashed). Their intersections 
represent periodic windows, band merging and laminar states. The vertical dotted lines show 
a choosing of points of band merging and laminar behaviour (a = 3.678, 3.727, 3.752, 3.791, 
3.877, 3.927). 



A useful tool for studying the chaotic behavior are the recursively formed super- 
track functions 



Si+i{a) = asi{a){l - Si{a)), so(a) = ^ 



(6) 



which represent the functional dependence of stable states |]33|1. The intersection of 
Si{a) with Si+j{a) indicates the occurrence of a j-period cycle and the intersection of 
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Si{a) with the fixed-point (1 — 1/a) of Eq. ^indicates the point of an unstable singular- 
ity, i. e. laminar behavior (Fig. |I|, intersection points are marked with dotted lines). 
For each a we compute a time series of the length N = 2, 000. In order to exclude 
transient responses we use the last 1, 000 values of these data series in the following 
analysis. 

We compute the RP after embedding the time series with a dimension of m = 1, 
a delay of r = 1 and a cut-off distance of e = 0.1 (in units of the standard deviation 
a). Since the considered example is a one-dimensional map, m = 1 is sufficient. In 
general, a too small embedding leads to false recurrences, which are expressed in 
numerous vertical structures and diagonals from the upper left to the lower right 
corner [ |30| 1. In contrast, an over-embedding should theoretically not distort the re- 
constructed phase trajectory. Whereas false recurrences and over-embedding do not 
strongly influence the measures based on diagonal structures [|30|1, the measures 
based on vertical structures are, in general, much more sensitive to the embedding. 
This is due to the fact, that the embedding method causes higher order correlations in 
the phase-space trajectory, which will be of course visible in the RP. A theoretical and 
more detailed explanation of this effect within the analysis of RPs is in preparation 
and beyond the scope of this article. For the logistic map, however, an increasing of m 
slightly amplifies the peaks of the vertical based complexity measures (up to m = 3), 
but it does not change the result significantly. The cut-off distance e is selected to be 
10 percent of the diameter of the reconstructed phase space. Smaller values would 
lead to a better distinction of small variations (e. g. the range before the accumulation 
point consists of small variations). However, the recurrence point density decreases 
in the same way and thus the statistics of continuous structures in the RP becomes 
soon insufficient. Larger values cause a higher recurrence point density, but a lower 
sensitivity to small variations. 



A. Recurrence Plots of the Logistic Map 

For various values of the control parameter a we obtain RPs, which already exhibit 
specific features (Fig. §). Periodic states (e. g. in the periodic window of length three 
at a = 3.830) cause continuous and periodic diagonal lines in the RP of a width of one. 



7 



There are no vertical or horizontal lines (Fig. |2|A). Band merging points and other 
cross points of supertrack functions (e.g. a = 3.720, Fig. gC) represent states with 
short laminar behavior and cause vertically and horizontally spread black areas in 
the RP. The band merging at a = 3.679 causes frequent laminar states and therefore 
a lot of vertically and horizontally spread black areas in the RP (Fig. 0B). Fully de- 
veloped chaotic states (a = 4) cause a rather homogeneous RP with numerous single 
points and rare short diagonal or vertical lines (Fig. §D). 

B. Complexity Measures of the Logistic Map 

Now we compute the known RQA measures A, L^ax and in addition (L) (average 
length of diagonal lines) and our measures A, Vmux and T for the entire RP of each 
control parameter a. As expected, the known RQA measures A, L^ax and (L) clearly 
detect the transitions from chaotic to periodic sequences and vice versa (Fig. |3|A, C, 
E) [^3p. However, it seems that one cannot get more information than the periodic- 
chaotic/ chaotic-periodic transitions. Near the supertrack crossing points (band merg- 
ing points included), e.g. a = 3.678, 3.791, 3.927, there are no significant indications in 
these RQA measures. They clearly identify the bifurcation points (periodic-chaotic/ 
chaotic-periodic transitions), without, however, finding the chaos-chaos transitions 
and the laminar states. 

Calculating the vertical based measures A and T, we are able to find the periodic- 
chaotic/ chaotic-periodic transitions and the laminar states (Fig. ^B, F). The occur- 
rence of vertical lines starts shortly before the band merging from two to one band at 
a = 3.678... 

For smaller a-values the consecutive points jump between the two bands and it is 
therefore impossible to obtain a laminar behavior. A longer persistence of states is not 
possible until all bands are merged. However, due to the finite range of neighborhood 
searching in the phase space, vertical lines occur before this point. 

Vertical lines occur much more frequently at supertrack crossing points (band 
merging points included), than in other chaotic regimes, which is revealed by A 
(cf Fig. ^B, again, supertrack crossing points are marked with dotted lines). As 
in the states before the merging from two to one band, vertical lines are not found 
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FIG. 2: Recurrence Plots (RP) of the logistic map for various control parameters a, near 
different qualitative changes: periodic-3-window a = 3.830 (A), band merging a = 3.679 (B), 
supertrack intersection a = 3.720 (C) and chaos (exterior crisis) a = 4 (D); with embedding 
dimension m = 1, time delay r = 1 and distance cutoff e = 0.1a. 
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FIG. 3: Selected RQA parameters A, Lmax and (L) and the new measures A, Vmax and T. 
The vertical dotted lines show some of the points of band merging and laminar behavior 
(cf Fig. |l|), whereby not all of them have been marked. Whereas A (A), L„iax (C) and (L) 
(E) show periodic-chaotic/ chaotic-periodic transitions (maxima), A (B), Vmax (D) and T (F) 
exhibit in addition to those transitions (minima) chaotic-chaotic transitions (maxima). The 
differences between A and Vmax are caused by the fact, that A measures only the amount of 
laminar states, whereas Vmax measures the maximal duration of the laminar states. Although 
some peaks of Vmax and T are not at the dotted lines, they correspond with laminar states 
(not all can be marked). 



within periodic windows, e.g. a = 3.848. The mean of the distribution of v is the in- 
troduced measure T (Fig. |3|F). It will vanish if a is smaller than the point of merging 
from two to one band. T increases at those points where more low ordered super- 
track functions are crossing (Fig. ^F). This corresponds to the occurrence of laminar 
states. Although Vmax also reveals laminar states, it is quite different from the other 
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two measures, because it gives the maximum of all of the durations of the laminar 
states. However, periodic states are also associated with vanishing T and Vmax- 

Hence, the vertical length based measures yield periodic-chaotic/ chaotic-periodic 
as well as chaos-chaos transitions (laminar states). 

We have also computed A, Vmax and T for the logistic map with transients using 
the same approach as described in | |23p . The qualitative statement of the measures is 
the same as above. 



V. APPLICATION TO HEART RATE VARIABILITY DATA 

Heart rate variability (HRV) typically shows a complex behavior and it is difficult 
to identify disease specific patterns | |34p . A fundamental challenge in cardiology is to 
find early signs of ventricular tachyarrhythmias (VT) in patients with an implanted 
cardioverter-defibrillator (ICD) based on HRV data ^, ^ Therefore stan- 
dard HRV parameters from time and frequency domain QS^], parameters from sym- 
bolic dynamics |I^ as well as the finite-time growth rates were applied to 
the data of a clinical pilot study [^]. Using two nonlinear approaches, we have re- 
cently found significant differences between control and VT time series based mainly 
on laminar phases in the data before a VT. Therefore the aim of this investigation is 
to test whether our RP approach is suitable to identify and quantify these laminar 
phases. 

The defibrillators used in the study cited (PCD 7220/7221, Medtronic) are able to 
store at least 1000 beat-to-beat intervals prior to the onset of VT (10 ms resolution), 
corresponding to approximately 9-15 minutes. We reanalyze these intervals from 17 
chronic heart failure ICD patients just before the onset of a VT and at a control time, 
i. e. without a following arrhythmic event. Time series including more than one non- 
sustained VT, with induced VT's, pacemaker activity or more than 10 % of ventricular 
premature beats are not considered in this analysis. Some patients had several VT's; 
we finally had 24 time series with a subsequent VT and the respective 24 control 
series without a life-threatening arrhythmia. In order to analyze only the dynamics 
occurring just before a VT, the beat-to-beat intervals of the VT itself at the end of the 
time series are removed from the tachograms. 
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We calculate all standard RQA parameters described in [ J24y as well as the new 
measures laminarity A, trapping time T and maximal vertical line length Vmax (in 
similarity to the maximal diagonal line length L^ax) for different embedding dimen- 
sions m and nearest neighbouring radii e. We find differences between both groups 
of data for several of the parameters mentioned above. However, the most significant 
parameters are Vmax and Lmax for rather large radii (Tab. The vertical line length 
Vmax is more powerful in discriminating both groups than the diagonal line length 
Lmax, as can be recognized by the higher p-values for Vmax (Tab. ||). Figure § gives 
a typical example of the recurrence plots before a VT and at a control time with an 
embedding of 6 and a radius of 110. The RP before a life-threatening arrhythmia is 
characterized by large black rectangles (Vmax = 242 here), whereas the RP from the 
control series shows only small rectangles (Vmax = 117). 




123456789 1 23458789 

Time [minj Time [min] 



FIG. 4: Recurrence plots of the heart beat interval time series at a control time (A) and 
before a VT (B) with an embedding of 6 and a radius of 110. The RP before a life-threatening 
arrhythmia is characterized by big black rectangles whereas the RP from the control series 
shows only small rectangles. 
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TABLE I: Results of maximal diagonal and vertical line length shortly before VT and at 
control time, nonparametric Mann- Whitney U-test, p - significance; * -p < 0.05; ** - p < 0.01; 
n. s. - not significant p > 0.05) 

m e VT Control p 



Maximal diagonal line length Lmax 

3 77 396.6±253.8 261.5±156.6 n. s. 

6 110 447.6±269.1 285.5±160.4 * 

9 150 504.6±265.9 311.6±157.2 * 

12 170 520.7±268.8 324.7±180.2 * 

Maximal vertical line length Vmax 

3 77 261.4±193.5 169.2±135.9 * 

6 110 283.7±190.4 179.5±134.1 ** 

9 150 342.4±193.6 216.1±137.1 ** 

12 170 353.5±221.4 215.1±138.6 ** 



VI. SUMMARY 

We have introduced three new recurrence plot (RP) based measures of complexity, 
the laminarity A, the trapping time T and the maximal length of vertical structures in 
the RP Vmax- These measures of complexity have been applied to the logistic map and 
heart rate variability data. In contrast to the known RQA measures ([^3[], [^]), which 



are able to detect transitions between chaotic and periodic states (and vice versa), our 
new measures enable us to identify laminar states too, i. e. chaos-chaos transitions. 
These measures are provided by the vertical lines in recurrence plots. The occurrence 
of vertical (and horizontal) structures is directly related to the occurrence of laminar 
states. 

The laminarity A enables us generally to detect laminar states in a dynamical 
system. The trapping time T contains information about the frequency of the laminar 
states and their lengths. The maximal length Vmax reveals information about the 
time duration of the laminar states thus making the investigation of intermittency 
possible. 
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If the embedding of the data is too small, it will lead to false recurrences, which is 
expressed in numerous vertical structures and diagonals perpendicular to the main 
diagonal. Whereas false recurrences do not influence the measures based on diagonal 
structures, the measures based on vertical structures are sensitive to it. 

The application of these measures to the logistic equation for a range of vari- 
ous control parameters has revealed points of laminar states without any additional 
knowledge about the characteristic parameters or dynamical behavior of the specific 
systems. Nevertheless, A, Vmax and T are different in their magnitudes. Further 
investigations are necessary to understand all relations between the magnitudes of 
Vmax and the recognized chaos-chaos transitions. 

The application of the new complexity measures to the ICD stored heart rate data 
before the onset of a life-threatening arrhythmia seems to be very successful for the 
detection of laminar phases thus making a prediction of such VT possible. The dif- 
ferences between the VT and the control series are more significant than in [|26[]. 
However, two limitations of this study are the relative small number of time series 
and the reduced statistical analysis (no subdivisions concerning age, sex and heart 
disease). For this reason, our results should to be validated on a larger data base. 
Furthermore, this investigation could be enhanced for tachograms including more 
than 10% ventricular premature beats. In conclusion, this study has demonstrated 
that the RQA based complexity measures could play an important role in the predic- 
tion of VT events even in short term HRV time series. 

Many biological data contain epochs of laminar states, which can be detected and 
quantified by the RP based measures. We have demonstrated differences between the 
measures based on the vertical and the diagonal structures and therefore we suggest 
the use of the method proposed in this article in addition to the traditional measures. 

A download of the Matlab implementation is available at: www.agnld.uni- 
potsdam.de/~marwan. 
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